Задача. Рассматриваются временные ряды NDVI. Эти данные обычно содержат шумовую компоненту, которая появляется вследствие влияния облачности, состояния атмосферы и т.п. Поэтому если построить график зависимости NDVI от времени, получится не красивая гладкая кривая, а сложная ломаная линия, в которой найти закономерность довольно сложно. Поэтому возникает задача очистки рядов NDVI от шумов.
Как читать. Начать стоит с раздела Теории.
Следующий за теорией раздел Реализация вначале содержит подраздел, в котором рассматривается вопрос о программной реализации фильтра, этот подраздел содержит много технических деталей, но может быть важен для понимания того, как функционирует описанная система. Последняя часть раздела "Реализация" содержит пошаговое выполнение алгоритма, сопровождаемое выводом получаемой на каждом шаге информации. Думаю, что для ознакомления с методом стоит прочитать этот раздел сразу после раздела теории, а технические детали оставить на потом.
Раздел Автоматизация содержит функции, готовые для использования в скриптах.
В конце находится раздел, посвященный исследованию поведения алгоритма на модельных данных и раздел, в котором обрабатываются реальные данные NDVI.
Методов обработки зашумленных данных NDVI много, здесь за основу процедуры взята статья: Chen, Jin, et al. "A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter." Remote sensing of Environment 91.3 (2004) 332-344..
В основе процедуры лежит фильтр Савицкого-Голая (Savitzky–Golay), применяемый для сглаживания изменяющихся во времени сигналов. Суть его состоит в том, что в скользящем окне шириной $2m+1$ исходная кривая аппроксимируется полиномом степени $d$ (например, параболой: $d=2$), коэффициенты которого расчитываются тем или иным способом (например, методом наименьших квадратов). Для центральной точки окна по этому полиному вычисляется значение, которое и будет результирующим значением отфильтрованной кривой в данной точке. Затем окно сдвигается на одну позицию и строится новая точка --- прогнозное значение кривой. Процедура повторяется до тех пор, пока не закончится исходный ряд.
Если же говорить об обработке серий NDVI, то у них есть особенности поведения, которые стоит учитывать. Одна из основных особенностей состоит в том, что шумовая компонента (возникающая вследствие облачности и других атмосферных явлений) на участках, покрытых растительностью, обычно понижает значение NDVI. Поэтому при фильтрации данных стоит находить кривую, огибающую график временной серии NDVI, по верхней границе, а не проходящую по центру серии.
Обозначения: дан временной ряд $(t_i, N_i)$, где $t_i$ --- дата $N_i$ --- значение NDVI в момент $t_i$. Дополнительно к $N_i$ у исследователя может быть ряд $F_i$ --- маска облачности в момент $t_i$.
Согласно статье обработка происходит по следующему алгоритму:
In [1]:
import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tools import get_ndvi_frame
from tools import get_modis_ndvi_frame
In [2]:
frame = get_ndvi_frame('data/Landsat/point_1_band3.csv', 'data/Landsat/point_1_band4.csv')
frame.head()
Out[2]:
Построим график NDVI
In [3]:
plt.plot(frame.DayNum, frame.NDVI)
Out[3]:
График перегружен информацией. Построим график за 1986-1986 год, чтобы рассмотреть некоторые детали.
In [4]:
selection = frame.ix[(frame.Year==1986)| (frame.Year==1987)]
print selection.head()
# Подготовим метки для показа на графике
# (отложим по оси абсцисс не все дни, иначе
# среди меток будет каша)
size = len(selection)
indexes = range(0, size, 2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(selection.DayNum, selection.NDVI)
plt.xticks(selection.DayNum[indexes], selection.Day[indexes], rotation='vertical')
plt.show()
Из этого графика более-менее видно, как ведет себя индекс в течении пары лет. В дальнейшем для примера будем работать с этой выборкой.
Построим фильтр на базе фильтра Савицкого-Голая. В scipy.signal реализован этот фильтр, но, к сожалению, в нем предполагается, что сигнал задан через равные интервалы (см. исходный код фильтра). Таким образом, если пользоваться готовым фильтром, то мы должны предварительно произвести интерполяцию NDVI так, чтобы временной ряд был равноотстоящим. Можно модифицировать код так, чтобы он мог работать с неполными (неравноотстоящими рядами). Проблема такого подхода в том, вычислительная сложность значительно возрастает.
Пока решил идти вторым путем. За основу функции взят исходный код http://wiki.scipy.org/Cookbook/SavitzkyGolay и модифицирован для работы неполными рядами.
Исходный код был такой (с небольшими упрощениями кода -- выкинул дифференцирование, которое нам не нужно):
In [5]:
def savitzky_golay(y, window_size, order):
r"""Smooth data with a Savitzky-Golay filter.
Parameters
----------
y : array_like, shape (N,)
the values of the time history of the signal.
window_size : int
the length of the window. Must be an odd integer number.
order : int
the order of the polynomial used in the filtering.
Must be less then `window_size` - 1.
Returns
-------
ys : ndarray, shape (N)
the smoothed signal.
"""
try:
window_size = np.abs(np.int(window_size))
order = np.abs(np.int(order))
except ValueError, msg:
raise ValueError("window_size and order have to be of type int")
if window_size % 2 != 1 or window_size < 1:
raise TypeError("window_size size must be a positive odd number")
if window_size < order + 2:
raise TypeError("window_size is too small for the polynomials order")
order_range = range(order+1)
half_window = (window_size -1) // 2
# precompute coefficients
b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
m = np.linalg.pinv(b).A[0]
# pad the signal at the extremes with
# values taken from the signal itself
firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
y = np.concatenate((firstvals, y, lastvals))
return np.convolve( m[::-1], y, mode='valid')
Создадим функцию, которая реализует процедуру для неполных рядов:
In [6]:
def savitzky_golay2(t, y, window_size, order):
r"""Smooth data with a Savitzky-Golay filter.
Parameters
----------
t : array_like, shape (N,)
the values of the time.
y : array_like, shape (N,)
the values of the time history of the signal.
window_size : int
the length of the window. Must be an odd integer number.
order : int
the order of the polynomial used in the filtering.
Must be less then `window_size` - 1.
Returns
-------
ys : ndarray, shape (N)
the smoothed signal.
"""
try:
window_size = np.abs(np.int(window_size))
order = np.abs(np.int(order))
except ValueError, msg:
raise ValueError("window_size and order have to be of type int")
if window_size % 2 != 1 or window_size < 1:
raise TypeError("window_size size must be a positive odd number")
if window_size < order + 2:
raise TypeError("window_size is too small for the polynomials order")
order_range = range(order+1)
half_window = (window_size -1) // 2
y = np.array(y) # It prevents caveats in pandas' indexing
# pad the signal at the extremes with
# values taken from the signal itself
firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
y = np.concatenate((firstvals, y, lastvals))
# Store the result
result = y.copy()
# Store distances between ticks
spaces = np.diff(t)
# copy first and last parst of 'spaces', to fit the lengths of 'spaces' and 'y' arrays
spaces = np.concatenate((spaces[:half_window], spaces, spaces[-half_window:]))
# precompute coefficients
for begin in range(len(y) - window_size):
# Running window
win = spaces[begin : begin + window_size].copy()
# Calculate the distances between the center point and the ticks:
win[half_window] = 0 # The center of the windows
first_half = win[:half_window].copy()
second_half = win[half_window:]
win[half_window:] = np.cumsum(second_half) # Distances to the right
win[:half_window] = -first_half[::-1].cumsum()[::-1] # Distances to the left
# Ordinary Least Squares:
b = np.mat([[k**i for i in order_range] for k in win])
m = np.linalg.pinv(b).A[0]
result[begin+half_window] = np.convolve( m[::-1], y[begin: begin + window_size], mode='valid')
return result[half_window: len(y) - half_window]
Проверим работу: создадим сигнал (с равными отсчетами) и прогоним его через обе функции. Для равных отсчетов результат должен быть одинаковым
In [7]:
order = 3
window_size = 7
t = np.arange(12)
y = t**2 + 10*np.random.randn(len(t))
print savitzky_golay(y, window_size, order)
print savitzky_golay2(t, y, window_size, order)
Данные по NDVI были импортированы в таблицу frame:
In [8]:
frame.head()
Out[8]:
Согласно используемой статье, необходимо произвести интерполяцию значений NDVI за те даты, которые приходились на облачную погоду.
Но во-первых, маски облачности у нас нет. А во-вторых, судя по всему, в статье описывался подход, когда фильтр работал с равными интервалами (а у нас есть фильтр savitzky_golay2, способный работать с неполными рядами).
Учитывая все сказанное, пропустим этот шаг.
Важно выбрать правильное значение $m$ --- ширины окна, в котором будем производить подгонку полинома под наши данные. Слишком маленькое значение приведет к тому, что в результирующей кривой останется излишинее количество деталей, но если выбрать большое $m$, то сглаживание будет избыточным. Авторы статьи рекомендуют начинать со значений $m$, равных 4--7. Аналогично со степенью полинома $d$. Попробуем начать с параболы (2-я степень) и постепенно увеличим до 4-й степени.
Таким образом "хорошая" пара $(m, d)$ подбирается экспериментально.
Будем анализировать не все данные, а только за пару лет. Произведем выборку и построим график NDVI.
In [9]:
selection = frame.ix[(frame.Year==1986)| (frame.Year==1987)].copy()
# Подготовим метки для показа на графике
# (отложим по оси абсцисс не все дни, иначе
# среди меток будет каша)
size = len(selection)
indexes = range(0, size, 2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(selection.DayNum, selection.NDVI, label='NDVI')
plt.xticks(selection.DayNum[indexes], selection.Day[indexes], rotation='vertical')
plt.legend(loc=2)
plt.show()
Прогоним через фильтр. Заметим, что в статье ширина окна измеряется в "половинках", а в функции savitzky_golay2 в "целых". Т.е. рекомендуемые авторами 4--7 единиц ширины означает 9-15 единиц в терминах функции savitzky_golay2. Попробуем начать с ширины равной 9-ти и порядка $d=2$
In [10]:
trend9 = savitzky_golay2(selection.DayNum, selection.NDVI, window_size=9, order=2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(selection.DayNum, selection.NDVI, label='NDVI')
plt.plot(selection.DayNum, trend9, label='Trend')
plt.legend(loc=2)
plt.xticks(selection.DayNum[indexes], selection.Day[indexes], rotation='vertical')
plt.show()
Возьмем окно в 15 отсчетов:
In [11]:
trend15 = savitzky_golay2(selection.DayNum, selection.NDVI, window_size=15, order=2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(selection.DayNum, selection.NDVI, label='NDVI')
plt.plot(selection.DayNum, trend15, label='Trend')
plt.legend(loc=2)
plt.xticks(selection.DayNum[indexes], selection.Day[indexes], rotation='vertical')
plt.show()
Лично мне в качестве первого прибилижения тренд, построенный на 15-ти отсчетах, нравится больше. Возьму за основу его.
In [12]:
trend = trend15
In [13]:
d = np.abs(selection.NDVI - trend)
dmax = np.max(d)
print 'dmax =', dmax
Для удобства дальнейшей работы добавим в таблицу selection еще две колонки: колонку, содержащую тренд, и колонку, содержащую веса.
In [14]:
selection['Trend'] = trend
print selection.head()
def get_w(row, dmax=dmax):
# Расчет весов для строки row
result = 1 if row['NDVI'] >= row['Trend'] else 1 - abs(row['NDVI']- row['Trend'])/dmax
return result
selection['W'] = selection.apply(get_w, axis=1)
print selection.head()
In [15]:
def get_n(row):
# Функция, в которой закодирована формула
result = row['NDVI'] if row['NDVI'] >= row['Trend'] else row['Trend']
return result
selection['N'] = selection.apply(get_n, axis=1)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(selection.DayNum, selection.NDVI, label='NDVI')
plt.plot(selection.DayNum, selection.N, label='New')
plt.legend(loc=2)
plt.xticks(selection.DayNum[indexes], selection.Day[indexes], rotation='vertical')
plt.show()
In [16]:
trend = savitzky_golay2(selection.DayNum, selection.N, window_size=13, order=2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(selection.DayNum, selection.NDVI, label='NDVI')
plt.plot(selection.DayNum, selection.Trend, label='Old Trend')
plt.plot(selection.DayNum, trend, label='New Trend')
plt.legend(loc=2)
plt.xticks(selection.DayNum[indexes], selection.Day[indexes], rotation='vertical')
plt.show()
selection.Trend = trend
Новый ряд действительно проходит выше, чем предыдущий.
In [17]:
selection['F'] = abs(selection['NDVI'] - selection['N'])*selection['W']
print selection.head()
f1 = selection.F.sum()
print 'F1 =', f1
Фукнкция фильтрации была определена выше.
Алгоритм, описанный в статье представлен в виде функции. На входе ожидается фрейм данных в формате, возвращаемом функцией импорта. Т.е. там жестко зашиты названия колонок (можно от этого избавиться, но, думаю, пока это не важно --- все равно итоговый модуль будет иным).
In [18]:
def fit_ndvi(frame, window_size=15, order=2):
'''Chen, Jin, et al. "A simple method for reconstructing a
high-quality NDVI time-series data set based on the
Savitzky–Golay filter."
Remote sensing of Environment 91.3 (2004) 332-344.
'''
def get_w(row, dmax):
# Расчет весов для строки row
result = 1 if row['NDVI'] >= row['Trend'] else 1 - abs(row['NDVI']- row['Trend'])/dmax
return result
def get_n(row):
# Функция, в которой закодирована формула
result = row['NDVI'] if row['NDVI'] >= row['N'] else row['N']
return result
# step 2
trend = savitzky_golay2(frame.DayNum, frame.NDVI, window_size=window_size, order=order)
frame['Trend'] = trend
# step 3
d = np.abs(frame.NDVI - trend)
dmax = np.max(d)
frame['W'] = frame.apply(get_w, axis=1, args=(dmax,))
iteration = 1
Old_F = 1000000 # A big number
Current_F = Old_F - 1
frame['N'] = frame['NDVI']
while Current_F < Old_F:
print 'iteration:', iteration, 'Current_F', Current_F
# step 4
frame['N'] = frame.apply(get_n, axis=1)
# step 5
# decrease window_size by t
t = min(8, 2 * iteration) # Just a test
trend = savitzky_golay2(frame.DayNum, frame.N, window_size=window_size-t, order=2)
frame.N = trend
frame['F'] = abs(frame['NDVI'] - frame['N'])*frame['W']
Old_F = Current_F
Current_F = frame.F.sum()
print 'New Current_F is assigned to', Current_F
iteration += 1
return frame
В этом разделе посмотрим, как ведет себя алгоритм на модельных данных. Для этого:
In [19]:
def get_signal(years=2):
T = years * 365 # Число дней, которые мы будем анализировать
bias = 40
day = np.arange(bias, T + bias, 10) # возьмем не все, а шагом в 10 дней
# Не важно, какая формула, главное чтобы было "похоже"
ndvi = np.maximum(np.sin(2*years*pi*day/T - bias) + 0.7, 0)
ndvi = (ndvi)**(0.3)
ndvi = ndvi/np.max(ndvi)
return pd.DataFrame(data = {'DayNum': day, 'NDVI': ndvi})
frame = get_signal()
plot(frame.DayNum, frame.NDVI)
Out[19]:
Создадим функцию, которая будет зашумлять сигнал и частично удалять данные из него. Функция принимает параметр drop_perc -- процент точек, которые будут удалены из сигнала, параметр noise_perc -- процент точек, которые будут зашумлены, и параметр noise_level --- значение, которому будет приравнена шумовая часть сигнала.
In [20]:
def get_noisy_signal(drop_perc=50, noise_perc=50, noise_level=0.3, years=2):
np.random.seed(0)
signal = get_signal(years=years)
size = len(signal)
all_index = signal.index.tolist()
noisy_count = int(size*noise_perc/100)
noisy_index = np.random.choice(all_index, size=noisy_count)
signal.NDVI[noisy_index] = noise_level + np.random.randn(noisy_count)/10
drop_count = int(size*drop_perc/100)
drop_index = np.random.choice(all_index, size=drop_count)
signal.NDVI[drop_index] = nan
return signal.dropna()
frame = get_noisy_signal()
plot(frame.DayNum, frame.NDVI)
Out[20]:
Прогоним через алгоритм слабозашумленный сигнал.
In [21]:
frame = get_noisy_signal(drop_perc=10, noise_perc=10)
plot(frame.DayNum, frame.NDVI)
Out[21]:
In [22]:
frame = get_noisy_signal(drop_perc=10, noise_perc=10)
filtered = fit_ndvi(frame)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(filtered.DayNum, filtered.NDVI, 'ro--', label='Noisy Points')
plt.plot(filtered.DayNum, filtered.N, label='Filtered')
plt.legend(loc=2)
plt.show()
Как видим, слабозашумленный сигнал воспроизводится достаточно точно. Однако в тех местах, где были горизонтальные участки сигнала, результирующая кривая излишне сгладила данные и эти участки искривились.
Возьмем сигнал более сложный:
In [23]:
frame = get_noisy_signal(drop_perc=25, noise_perc=25)
filtered = fit_ndvi(frame)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(filtered.DayNum, filtered.NDVI, 'ro--', label='Noisy Points')
plt.plot(filtered.DayNum, filtered.N, label='Filtered')
plt.legend(loc=2)
plt.show()
Здесь "горбы" исходной кривой восстановлены довольно точно. Но в промежутке между горбами отфильтрованная кривая задралась вверх. Частично это явление мы видили еще на предыдушем рисунке, а частично на него повлияло то, что на этом участке появилось несколько шумовых точек, повышающих уровень исходного сигнала. Поскольку в алгоритм заложена процедура "подтягивания" сглаженного сигнала к точкам наибольшего уровня, то видим, что в итоге промежуток между горбами поднялся с нуля до 0.35--0.4 единиц.
Усилим зашумление еще больше:
In [24]:
frame = get_noisy_signal(drop_perc=50, noise_perc=50)
filtered = fit_ndvi(frame)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(filtered.DayNum, filtered.NDVI, 'ro--', label='Noisy Points')
plt.plot(filtered.DayNum, filtered.N, label='Filtered')
plt.legend(loc=2)
plt.show()
На этом рисунке ничего принципиально нового не появляется: видим, что горбы восстановлены с достаточной точностью (хотя погрешность, конечно, возрасла), а промежуток между ними оказался излишне задранным вверх (в большей степени, чем на предыдущем рисунке).
И, наконец, чтобы поиздеваться вволю над машиной, подадим сплошной шум:
In [25]:
frame = get_noisy_signal(drop_perc=70, noise_perc=70)
filtered = fit_ndvi(frame)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(filtered.DayNum, filtered.NDVI, 'ro--', label='Noisy Points')
plt.plot(filtered.DayNum, filtered.N, label='Filtered')
plt.legend(loc=2)
plt.show()
Попробуем посмотреть на "многолетние" циклы:
In [26]:
frame = get_noisy_signal(drop_perc=10, noise_perc=10, years=8)
filtered = fit_ndvi(frame)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(filtered.DayNum, filtered.NDVI, 'ro--', label='Noisy Points')
plt.plot(filtered.DayNum, filtered.N, label='Filtered')
plt.legend(loc=2)
plt.show()
Слабо зашумленный сигнал восстанавливается вполне приемлемо. Проблемы все те же, что и замеченные ранее, новых не добавилось.
Добавим шумов:
In [27]:
frame = get_noisy_signal(drop_perc=50, noise_perc=50, years=8)
filtered = fit_ndvi(frame)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(filtered.DayNum, filtered.NDVI, 'ro--', label='Noisy Points')
plt.plot(filtered.DayNum, filtered.N, label='Filtered')
plt.legend(loc=2)
plt.show()
В целом восстановление сигнала на модельных данных идет удовлетворительно.
Основная замеченная проблема: завышение уровня сигнала на участках с низким NDVI. В случае низкого NDVI (зимнее время, участки не покрытые растительностью) весьма вероятно появления точек со значениями, превышающими "нормальный" уровень NDVI. В этих случаях процедура будет завышать показания.
Уменьшить эту проблему можно, если каким-либо образом оценить, является ли данная точка шумовой или нет. Видимо, один из самых простых способов --- использовать маску облачности (как и было рекомендовано в исходной статье).
In [28]:
frame = get_ndvi_frame('data/Landsat/point_1_band3.csv', 'data/Landsat/point_1_band4.csv')
selection = frame.ix[(frame.Year==1985) | (frame.Year==1986)| (frame.Year==1987) | (frame.Year==1988)].copy()
# Подготовим метки для показа на графике
size = len(selection)
indexes = range(0, size, 2)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(selection.DayNum, selection.NDVI)
plt.xticks(selection.DayNum[indexes], selection.Day[indexes], rotation='vertical')
plt.show()
In [29]:
filtered = fit_ndvi(selection)
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(filtered.DayNum, filtered.NDVI, 'ro--', label='Noisy Points')
plt.plot(filtered.DayNum, filtered.N, label='Filtered')
plt.legend(loc=2)
plt.show()
Основные проблемы алгоритма, которые были замечены на модельных данных никуда не изчезли. Алгоритм все так же завышает значения в зимний период. В тех случаях, когда зимних данных мало, и алгоритму по сути нечего завышать, мы видим почти линейную интерполяцию. Если данные за зимний период есть, то они результат пройдет по верхней (в данном случае шумовой) границе значений.
В качестве входных данных возьмем MOD13 --- 16-ти дневный продукт. Этот продукт содержит информацию о качестве полученного значения NDVI, которую можно учитывать в алгоритме фильтрации при назначении весов. Но пока будем использовать алгоритм "как есть", без адаптации.
In [30]:
datafile = 'data/MODIS/point_5_NDVI.csv'
frame = get_modis_ndvi_frame(datafile)
plot(frame.DayNum, frame.NDVI)
Out[30]:
In [31]:
filtered = fit_ndvi(frame.copy())
fig = plt.figure(figsize=(12,6), dpi=100)
plt.plot(filtered.DayNum, filtered.NDVI, 'ro--', label='Noisy Points')
plt.plot(filtered.DayNum, filtered.N, label='Filtered')
plt.legend(loc=2)
plt.show()
За счет более гладких (и более равномерно распределенных) входных данных MODIS алгоритм фильтрации работает значительно лучше, чем на Landsat. Тем не менее, на графиках отчетливо видно, что этот алгоритм "удлинняет" вегетативный период -- в осеннее и весеннее время отфильтрованные значения изменяются менее резко, чем реальные данные. Летний же период восстанавливается относительно нормально.